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We introduce simple prescriptions of the Yukawa potential to describe the effect of size polydis- 
persity and macroion shielding effect in charged colloidal systems. The solid-liquid phase boundaries 
were presented with the Lindemann criterion based on molecular dynamics simulations. Compared 
with the Robbins-Kremer-Grest simulation results, a deviation of melting line is observed at small 
A, which means large macroion screening length. This deviation of phase boundary is qualitatively 
consistent with the simulation result of the nonlinear Poisson-Boltzmann equation with full many- 
body interactions. It is found that this deviation of the solid-liquid phase behaviour is sensitive to 
the screening parameter. 

PACS: 82.70.Dd, 82.70.-y, 05.20.Jj. 

Charged colloidal suspensions have received much attention in recent years because of their abundant structural 
and rheological behaviour. As a classical complex many-body system, they are also an impressive model system since 
it is easy to be verified by various optical techniques. M The classical explicit expression of effective pair interaction of 
one-component model of colloidal systems is the Dejaguin-Landu-Verwey-Overbeek (DLVO) potential, which means 
the interaction between an isolated pair of spherical macroions with finite size. The DLVO potential consists of 
an electrostatic and the van der Waals part. I 2 ' In general, it is able to manage with Yukawa or screened Coulomb 
form which only retains the electrostatic part. Furthermore, the phase diagram of Yukawa systems was presented by 
computer simulations. t 3 ' 4 ' The phase transition line from order (bcc or fee) to disorder (liquid) is approximately given 
by a straight line, Tm = 0.00246 + 0.000274 • A, where Tm is the rescaled melting temperature and A is a measure of 
verse rescaled screening parameter. 

However, in actual colloidal suspensions, the particles are not identical in size, or they exist size polydispersity. In 
recent experiments Yamanaka et al. studied the reentrant crystalline-liquid transition induced by surface charge. ' 5_7 ' 
A visible discrepancy of phase boundary between size polydispersities 2% and 4% of latex is observed. On the 
other hand, Brunner et al. experimentally investigated the effective interparticle potential in two-dimensional charge- 
stabilized colloidal system at various particle densities. I 8-10 ! It is found that the effective pair colloid-colloid potential 
is in agreement with the Yukawa form of the DLVO theory only at low density, but not at high density. A suggestion, 
the so-called macroion shielding, has been put forward to explain the density-dependent pair interaction. 

It is well-known that the macroscopic properties of charged colloidal suspension are determined by microscopic 
interaction among the constituents of this system. The main purpose of the present work is to investigate the influence 
of macroscopical solid-liquid phase behaviour caused by size polydispersity and macroion shielding. Through picturing 
these effects as simple prescriptions of Yukawa pair interaction, we present the melting lines based on molecular 
dynamics (MD) simulation. A very similar deviation of the solid-liquid phase boundary is observed at large screening 
length of macroion. This indicates that the screening parameter plays an important role in the macroscopic properties 
of charge-stabilized colloidal suspensions. 

In the DLVO description of charged colloidal system of N macroions with number density p, the resulting effective 
pair potential between macroions is of the Yukawa shape with length scale a = p~ 1//3 , 

U(r)=A eM -* r/a \ (1) 
r/a 

where r is the center-center separation between two particles; A = na with k being the Debye-Hiikel screening 
parameter; and A — (Z*e) 2 / (ea) is the prefactor to describe the interaction strength with Z* being the charge of 
particle, e is the charge of an electron in electrostatic units, and e is the background solvent dielectric constant. It is 
noted that in this description of prefactor, we have incorporated the size correction factor directly into the effective 
charge. Clearly, the Yukawa picture of colloidal systems only depends on the rescaled separation between particles 
r/a and the parameters A and A. 
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In the course of our work on melting lines, we follow the simulation scheme of Robbinset al. I 3 ' 4 ! It is worthwhile 
to normalize the parameters with typical phonon energy. As a general rule, the rescaled temperature is defined by 

f = K B T/ (ma 2 u; 2 E ) , (2) 

where KbT is the thermal energy, m is the mass of particle, and uje is the Einstein frequency. Moreover, from the 
calculation of the total energy per particle at zero temperature (all particles are posited at the lattice sites), it is easy 
to characterize the prefactor of pair potential by Einstein frequency loe- 

In the sense of intuitive physics, the direct solution of studying the effect of size polydispersity is to deduce the 
fluctuation distribution of prefactor in the Yukawa potential from the size distribution of particles. However, it is 
a puzzle how to establish the dependence correlation between the distribution of effective surface charge and the 
size distribution. It is concluded that the different presuppositions, constant charge density and constant electrical 
potential on the surface of particle would lead to conflicting results. On the assumption of constant surface electrical 
potential, one could bring on the DLVO Yukawa pair potential, whether by analytic deduction from the mean spherical 
approximation of primitive model or by numerical simulation from complete Poisson-Boltzmann equation. I 11-14 ! 
On the other hand, for the assumption of constant surface charge density, it is proven from the hypernetted-chain 
approximation that the colloid-colloid interaction is not a simple Yukawa shape with neglect of van der Waals forces 
and it is possible to perform the attractive interaction at small separation under some conditions. I 15 ' 

In the viewpoint of theoretic works, with discarding the exact distribution function of charge numbers of particles 
caused by size polydispersity in charged colloidal systems, we have just introduced a simple prescription into the 
Yukawa potential, i.e., 

jj = A / cx P(~Ar/a) 
r/a 

A' = ,4(7+1), (3) 

where A is the prefactor in Eq. (1) with its value determined from Table 1 in Ref. [4] and 7 is a Gaussian random 
number with mean value zero and variance a . In Eq. (3), the effect of size polydispersity is approached by a Gaussian 
fluctuation of prefactor. It is not a accurate detailed description of size effect. Even so, this coarse-grained prescription 
is enough to extract very useful information about the effect of size polydispersity on the solid-liquid phase behaviour 
of charged colloidal suspensions. 

To simulate the solid-liquid phase boundary, in strictly speaking, one should calculate the free energies for solid and 
liquid states. In the past decades, there are many computer simulations for quantitative results and relative theoretic 
techniques about criterion of this phase transition, melting or freezing. All of these progresses are contributed by the 
advanced experiments and well-characterized colloidal systems. In this work, we are only interested in the solid-liquid 
phase transition, or melting line. Normally, one should be considered the Lindcmann criterion or the Hansen- Verlet 
rule in the calculation. The Lindemann criterion states that the ratio of the root-mean-square (rms) displacement 
\J (Su 2 ) and the characteristic interparticle distance or a of the solid at the Tm, the melting temperature, has a 
universal value. I 16 ! The Hansen- Verlet rule means that the first maximum of liquid structure factor has an amplitude 
of about 2.85 along the phase transition line. I 17 ' Both criteria would lead to the identical melting temperature in 
this simulation scheme. However, the rms displacement is easier to be calculated and gives the more precise values 
of Tm- ' 4 ' C onseq uently, the rms displacement is calculated to achieve the melting line according to the Lindemann 
criterion, yj {5u 2 )/a — 0.19, in all of our simulations. 

A very important simulation technique, constant temperature MD method, is employed to realize calculation of the 
rms displacement at expectant temperature. The so-called stochastic method that the experiment system is coupled 
with a head bath is used in our simulations. For every particle with mass m interacted by a total Yukawa potential 
U t from other particles, one can consider the equations of motion for Cartesian coordinate, I 4 ' 18 ' 19 ! 

d 2 v l dU t dr, 

where a friction force with coefficient 77 and a random force Ri (t) are contacted. The amplitude of random force 
is governed by the fluctuation dissipation theorem, (Rj (t\) Rj (£2)} = ^KsTriSijd (ti —t%), where T denotes the 
temperature of heat bath. This solution of constant MD simulation is easy not only to put into practice, but also to 
obtain excellent performance. In our simulations, the average relative error between the real stable temperature and 
the expectant temperature is less than 1%. 

All of our simulations was carried out with TV = 432 and 500 particles, for initial configurations of bec and fee 
crystals in cubic box with periodic boundary, respectively. The time step in simulations is At = O.OIte where 
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te = 2tt/loe is from the definition of energy scale. The time to reach equilibrium status for experimental systems is 
not longer than 10000 time steps, and the all runs is 100000 time steps. 

In Fig.l, we presented the melting lines of charged colloidal system with our simple prescription Eq. (3) for varied 
distributions. The open squares, circles and triangles are the melting points from our MD simulations according to 
the Lingdemann criterion for the variance of fluctuation a = 0.10, 0.05 and 0.01, respectively. The solid line is the 
RKG melting line and the others are the fitting curves. Obviously, at small A, or large Debye screening length, the 
deviation of solid-liquid phase boundary is observed. Moreover, the larger a from the size polydispersity becomes, 
the larger deviation of melting line becomes. For the larger values of A, however, it is just a very little shift of phase 
boundary compared with the RKG line. 

In a general way, the interaction between particles is characterized with the new length scale a, U a = t/oexp(— A). W 
In the experiment on phase diagram of charged colloidal suspensions, ^ it proves very good agreement between the 
experimental results and the RKG transition lines after using renormalized charge. For example, from Fig. 8(b) 
in Ref. [20], KBT/U a is about 0.15 for A = 5, but from the RKG prediction without renormalization, this value is 
0.1954. In our simulation result as shown in Fig.l, K B T/U a is 0.1937, 0.1883, and 0.1637 for a = 0.01, 0.05, and 0.10, 
respectively. 




FIG. 1. Melting lines of Yukawa systems with various fluctuations of prefactor. The solid line is the RKG melting line 
without fluctuation of prefactor. The length of cutoff for pair interaction r c = 3.07a for all phase boundaries. The squares, 
circles, and triangles indicate the values of Tm where the rms displacement reached 0.19 of a from our MD results. 



A qualitative rational explanation follows from the deviation in Fig.l at small A. The deviation of the solid-liquid 
phase boundary arises from the correlated fluctuations of the ion clouds around macroion. In a general case, this effect 
is not strong enough to influence the phase transition especially in monovalent electrolytes. However, this correlated 
effect has been enhanced by size polydispersity. As a result, the noticeable deviation is found at large polydispersity, 
such as c — 0.05,0.10 in Fig.l. For a = 0.01, the melting line is very close to the RKG line since the deviation at 
small A can be accepted after considering the simulation error. On the part of large A, the deviation has reduced to 
a very small shift because the ionic distribution about the central macroion decreases dramatically with distance and 
the polydispersity scarcely affects the correlated effect. As to the shift of phase boundary in Ref. [20], is is caused by 
the nonlinear effect in expansion of Possion equation. Anyway, in the physical nature, the all shifts of phase boundary 
is from the microscopic interaction although one root in the correlated fluctuation of ionic cloud and another comes 
from the nonlinear effect in expansion of Poisson equation. The detailed microscopic origins will be studied in the 
future works. 

All the above calculations are completed with the cutoff of pair interaction r c — 3.07a, as the same as for the 
RKG line, since the longer range contribution to the energy is close to an independent structure. M Acting up to 
the physical instinct, someone may image that the length of cutoff is longer, and the simulation result is better. 
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However, one should note that the Yukawa potential is only interaction between two particles with finite size, and is 
deduced with the presupposition of dilute limit of particle density. In fact, the three-body interaction and more order 
effects are obliterated in the calculation with this form of pair interaction. Based on the recent experimental results 
about effective pair interactions at different densities, I 8,9 ' we have proposed a simple prescription to the Yukawa 
potential, [21 " 231 

__. . . exp(— Xr/a) 
U(r) = A — ^— — — r <r c 
r I a 

= r > r c (5) 

with prefactor A : and a density-dependent cutoff r c oc a — p^ 1 ^ 3 to include the macroion shielding effect. 

Based on this simple modified pair potential, the solid-liquid phase boundary for r c — 1.77a is plotted in Fig. 2 as a 
dashed line and the open circles are the rms displacement calculated from our MD simulations with the Lindcrmann 
criterion. The solid line, f M = 0.00246 + 0.000274A, is the RKG melting line produced for r c = 3.07a. M It is easy 
to observe that a visible discrepancy is observed for small A, just like that in Fig.l. However, for a larger value of 
A, the melting points are almost identical to the RKG results and there are no shifts for all A. As interpretation 
mentioned above, the density-dependent truncation of pair interaction is resulted from the many-body effect in 
concentrated charge colloidal suspensions. It follows qualitatively that the macroscopical phase behaviour is only 
influenced noticeably for small A (or large Debye screening length) from Fig. 2. This result is consistent with the 
conclusion from simulations with a continuous Poisson-Boltzmann description which includes fully many-body effects 
in colloidal systems. I 21-23 ! 

To compare the above two simulation results, Fig.l and Fig. 2, we have observed a similar deviation of solid-liquid 
phase behaviour in the area of large screening length. However, in Fig.l, this discrepancy is from the fluctuation of 
interactions between particles, which is enhanced by the size polydispersity of colloidal particles. Figure 2 shows a 
result from the macroion shielding effect which is due to the many-body effect at high density. 




FIG. 2. Melting line of Yukawa systems with length of cutoff r c = 1.77a. There is no effect of size polydispersity added to 
experiment systems. The open circles are the values of Tm from the Lindemann criterion in our MD simulations. 



It is perhaps worthwhile to comment many-body attraction. It is well known that an attractive component in pair 
interaction can appear due to the correlated effect. As a straightforward result of strong correlation effect from large 
size polydispersity, the pair potential should include an attraction part. In our simulation for Fig.l, avoiding the 
difficulty to obtain the exact uniform effective pair interaction, the correlated fluctuation of ionic cloud is directly 
and approximately reflected by Gaussian distribution of the prefactor in the Yukawa pair potential. On the other 
hand, the charge-like attraction is obtained and it becomes visible with the increasing particle density. I 8,9 ! It is noted 
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that the effective part in pair potential is similar to the same part of the Yukawa potential after considering the 
density-dependent truncation of interaction length to include many-body effect. In a sense, we have conjectured that 
the deviation of macroscopic phase behaviour corresponds to appearance of many-body attraction especially at low 
salt concentration for high particle density or large size polydispersity. 

In conclusion, considering size polydispersity and macroion shielding effect as simple prescriptions of the Yukawa 
potential between colloidal particles, we have performed MD simulation results of solid-liquid phase behaviour. A 
visible deviation of phase boundary is observed at small A, or large screening length of macroions. It is found that this 
discrepancy is sensitive to the Debye screening length. Again, it is emphasized that the screening parameter plays a 
major role in the interpretation of macroscopic properties of charged colloidal systems. 

We gratefully acknowledge Jupei Yamanaka for useful and stimulating communications. 
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